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ABSTRACT 

Isolated regions of higher density populate the ISM on all scales - from molecular clouds, to 
the star-forming regions known as cores, to heterogeneous ejecta found near planetary nebulae 
and supernova remnants. These clumps interact with winds and shocks from nearby energetic 
sources. Understanding the interactions of shocked clumps is vital to our understanding of the 
composition, morphology, and evolution of the ISM. The evolution of shocked clumps are well- 
understood in the limiting "adiabatic" case where physical processes such as self-gravity, heat 
conduction, radiative cooling, and magnetic fields are ignored. In this paper we address the issue 
of evolution and convergence when one of these processes - radiative cooling - is included. 

Numeric convergence studies demonstrate that the evolution of an adiabatic clump is well- 
captured by roughly 100 cells per clump radius. The presence of radiative cooling, however, 
imposes limits on the problem due to the removal of thermal energy. Numerical studies which 
include radiative cooling typically adopt the 100-200 cells per clump radius resolution. In this 
paper we present the results of a convergence study for radiatively cooling clumps undertaken 
over a broad range of resolutions, from 12 to 1,536 cells per clump radius, employing adaptive 
mesh refinement (AMR) in a 2D axisymmetric geometry (2.5D). We also provide a fully 3D 
simulation, at 192 cells per clump radius, which supports our 2.5D results. We find no appreciable 
self-convergence at ^100 cells per clump radius as small-scale differences owing to increasingly 
resolving the cooling length have global effects. We therefore conclude that self-convergence is an 
insufficient criterion to apply on its own when addressing the question of sufficient resolution for 
radiatively cooled shocked clump simulations. We suggest the adoption of alternate criteria to 
support a statement of sufficient resolution, such as the demonstration of adequate resolution of 
the cooling layers behind shocks. We discuss an associated refinement criteria for AMR codes. 

Subject headings: hydrodynamics - ISMxlouds - shock waves 

1. Introduction 

The distribution of matter is not uniform across a wide range of scales in the universe. Within our own 
galaxy the ISM is not smooth and, in particular, density inhomogenities are found in molecular clouds, and 
within these clouds matter is distributed unevenly all the way down to the star-forming regions known as 
molecular cloud cores. Clumps of material exist on smaller scales as well. This heterogeneous distribution 
of matter is required, of course, for star and planet formation. On the other hand, energetic sources such as 
young stellar objects (YSOs), planetary nebulae (PNe), and supernovae inject kinetic energy back into their 
environments in the form of winds, jets, and shocks. It becomes a matter of central interest, therefore, to 
understand how clumps and winds, jets, and shocks interact. 
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One would expect a wind or shock coming from an astrophysical source to have some curvature. It 
follows that a shock sweeping over a spherical clump will have a radius of curvature which may or may not 
be comparable to that of the clump. For analytic simplicity many studies consider the "small clump" limit, 
in which the curvature of the shock is so large compared to the clump that it may be considered planar (see 



e.g. Klein et al. 1994 hereafter KMC94). Studies where this limit is not applied include those addressing 



questions such as the evolution of a heterogeneous envelope surrounding PNe (|Steffen fc Lopez 20041, or 



energetic jets impinging on a single clump (e.g. Raga et al. 2002; Fragile et al. |2004 ) or a host of clumps 
QYirak et al.]|2008 l. 



Initial studies of single clump/shock interaction focused on the early stages of the interaction, where 
the solution remained amendable to linear approximations. The evolution late in time, when the behavior 
becomes highly nonlinear, remains intractable from a purely analytic standpoint and therefore has benefitted 
greatly from numerical investigation — a review of the pioneering literature may be found in KMC94 or 
jludnenko et al. (20021. Illustrating the maturity of the field, a variety of physical processes have been 
included in the studies. KMC94 discussed systematically the evolution of a single, adiabatic, nonmagne- 
tized, non-thermally conducting shocked clump overrun by a planar shock in axisymmetry ("2.5D"). Similar 



simulations were carried out in three dimensions (3D) by Stone & Norman (19921. Magnetic fields can play 
an important part in the clump evolution — the role of field orientation was first investigated by | Jones | 
et al. (1996) and more recently by van Loo et al. (20071; Shin et al. (2008[). Smooth cloud boundaries (e.g. 



Nakamura et al. 20061, and systems of clumps (e.g. Poludnenko et al. |2002[ ) also have been studied. The 



similar problem involving clump-clump collisions has also received attention (Miniati et al. 1997 Klein & 



Woods 1998). Studies predominantly use an Eulerian mesh with a single- or two-fluid method to solve the 
inviscid Euler equations. One notable exception is Pittard et al. (20091, who use a "k — e" model to explicitly 



model the effects of sub-grid-scale turbulent viscosity by adding turbulence-specific viscosity and diffusion 
terms to the fluid equations. 



The role of radiative cooling has also been studied (e.g. Mellema et al.|2002 Fragile et al. 2004 Orlando 



et al. 20051. The energy loss associated with optically thin radiative cooling significantly modifies the evo- 



lution of the shocked clump. In the adiabatic regime, the shocked clump collapses initially before expanding 
in the direction perpendicular to the shock due to the increased thermal energy. For strong cooling, |Mellema| 



et al. (20021; Fragile et al. (20041 find that the clump not only does not reexpand but in fact fragments 



and collapses down to grid scales. Orlando et al. (2005), who consider the role of thermal conduction as 



well as magnetic fields and radiative losses, find that even for weaker cooling the evolution differs from the 
adiabatic case. There exists a broad range of velocities, temperatures, and densities that admit radiative 
cooling. Radiative cooling, therefore, is a significant effect which requires careful investigation. 

As with any numerical investigation the appropriate choice of spatial resolution is of critical importance. 
KMC94 suggested that ^100 cells per clump radius was required to adequately model the interaction. 
Authors subsequent to KMC94 typically have adopted resolutions in the range of 100-200 cells per clump 
radius when performing simulations in 2- or 2.5D. Three-dimensional simulations being more computationally 



expensive, authors in some cases choose lower resolution. A small 3D resolution study in Agertz et al. (20071, 
for example, investigated the behavior of a 3D adiabatic shocked clump with ^6, ~12, and ^25 cells per 
clump radius. In contrast, Shin et al. (2008) employed 120 cells per clump radius to model the effects of 



imposed magnetic fields on a single adiabatic shocked clump in 3D. Orlando et al. (2005) also addressed 



resolution by considering two 3D simulations, at 105 and 132 cells per clump radius, respectively, and found 
good agreement between the two. 



Is the choice of 100-200 cells per clump radius always appropriate in clump simulations when different 
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physical processes are included? Klein & Woods ( 1998) investigated clump-clump collisions in order to better 
understand the nonlinear thin shell instability (NTSI). The NTSI can result at the interface between two 
opposing flows, causing any perturbation at the interface to grow. When modelling adiabatic or radiatively 
cooling clumps, they found 128 cells per clump radius to be sufficient. When modelling isothermal clumps, 
however, they found that at this resolution the NTSI was damped, and resolutions 4-5 times higher were 
required to model the instability. The cause was the vanishing of the thin interaction layer behind isothermal 
shocks, which was reduced down to at most a few cells. This, combined with the specifics of their numerical 



scheme, resulted in an artificial numerical damping of the instability. Thus Klein & Woods (1998) were 
required to increase the resolution to counteract this effect. Such a requirement of adequate resolution 
applies not only to clump-clump collisions, of course. In simulations of environments where self-gravity plays 
a role, for example, typically the resolution requirements are described in terms of adequate resolution of 
the local Jeans length (see e.g. Truelove et al.fl998 ). It is always important that the relevant physics on the 



scales of interest be adequately resolved. 

This resolution/scale rule holds true for radiative cooling. The length scale of interest for radiative 
cooling behind shocks is the cooling length: the distance behind the radiative shock where the shocked 
material cools according to some cooling function oc n 2 A(T). As will be demonstrated, such a scale may 
not be adequately resolved in all cases previously considered. On the one hand, when not fully resolved, the 
cooling length will be artifically increased, implying that material a given distance behind a shock is at a 
higher temperature than it should be. On the other hand, if the cooling length is highly underresolved, the 



cooling will take place on of order the cell size, and the simulation will be essentially isothermal (Dyson & 



Williams 1997). If, as is the case in shocked clump interactions, multiple shocks interact with each other, 



this underresolution may lead to global incorrect results. 

In contrast to studies of adiabatic clumps, there has not been a convergence study undertaken for cooling 
clumps. A possible reason for this is that different cooling tables or multiphysics are employed in different 
studies, so the fine details may not be generalizeable. However, it remains worthwhile to study questions 
pertaining to cooling length as a necessary criterion for the resolution of radiatively cooling shocked clump 
simulations. 

We therefore seek to understand the effect of varying resolution in radiatively cooling shocked clump 
simulations. In particular, we address a sequence of linked questions: 1) how do the simulations vary, 
qualitatively and quantitatively, with resolution; 2) do the simulations converge in the same manner as 
adiabatic clump models; 3) since these simulations are performed in 2.5D, how might the conclusions extend 
to fully 3D simulations; and 4) what conclusions may be drawn in general about simulations of radiatively 
cooled clumps? 

In §[2] we describe the problem and physical scales of interest, and § [3] discusses our numerical methods. 
§ [4] highlights the results from the investigation, which are discussed in more quantitative detail in § [5] We 
conclude in § [6j 



2. Description of the Problem 

2.1. Evolution of a shocked clump 

The evolution of an adiabatic shocked clump progresses through several stages. The external shock will 
sweep across or "pass over" the clump radius in a time t sp = r c /v Sl where r c is the radius of the clump and 
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v s the shock velocity. When the shock hits the clump, a new shock wave is transmitted into the clump. This 
transmitted shock crosses the clump in a time known as the cloud- crushing time t cc (KMC94), 



X 1/2 r c 



(1) 



where v SjC is the transmitted (internal) shock velocity and x — Pc/ Pa, the clump-to- ambient density ratio. 
The approximation u s c ~ Vg/x 1 ^ 2 comes from considering approximate ram pressure balance for high-Mach 
number shocks, where p a v 2 — p c v 2 c - Since x > 1, then t sp < t CCl and the cloud will be enclosed in the hot 
postshock flow before the transmitted shock crosses the clump. No longer in pressure balance, the clump will 
begin to be compressed on all sides via new "transmitted shocks." When the transmitted shocks meet and 
interact, a rarefaction is produced which causes the clump to reexpand. Thereafter, the clump is destroyed 
and mixes with the postshock flow in a few cloud-crushing times. See KMC94 for a more detailed description 
on the evolution of an adiabatic shocked clump. How radiative cooling modifies this evolution is discussed 
in §|pl 



The two main sources of disruption of the clump at the clump/postshock interface are the Kelvin- 
Helmholtz (KH) and Ray lcigh- Taylor (RT) instabilities. The KH instability will occur due to the relative 
velocity between the postshock flow and the clump material. The growth rate of the KH instability in the 
case of x ^ 1 is £kh = X 1 ^ 2 / '(kv re i) ( Chandrasekhar 1961), where v re i is the relative velocity between the 
clump and postshock material. If we assume the relative velocity is nearly equal to the external postshock 
velocity (v ps = 3v s /4 for M»l), then the KH growth time is 



^kh 



4x 



1/2 



3fcv, 



ice 
kr r . 



(2) 



While the smallest wavelength instabilities will grow the fastest, it has been seen that the most disruptive 
modes are those of order the clump radius (KMC94, |Poludnenko et aL||2002 1 . Therefore, for the dominant 
modes for an adiabatic shock clump interaction will occur on timescales of £kh ~ t cc . 

The shocked clump also will be susceptible to the RT instability, as a lighter material (the postshock 
flow) is accelerating a heavier one (the clump). In the plane-parallel case, one would expect RT "bubbles" 
to move into the high-density clump and RT "spikes" to move into the lower-density postshock flow. For 
a spherical clump, the curvature of the clump leads to nonaxial flow which prevents the formation of the 
RT spikes. The RT growth time is tux = (gk) -1 ^ 2 (Chandrasekhar 19611, where g is the acceleration. As 
discussed in e.g. Fragile et al. (20041, we may approximate the acceleration as g ~ v s /t acc ~ v s /(x 1 ^ 2 t cc ), 
where t acc is the acceleration timescale. The acceleration timescale is approximated by considering the 
momentum transfer of postshock material onto the clump. With this approximation, the RT growth time is 



tRT 



tr. 



(3) 



(Ar c )V2 ' 

As with the KH instability, wavelengths of order the clump radius are seen to be the most disruptive. 
Therefore, as with the KH instability, for the dominant modes tux ~ t cc . By comparing the two instabilities 
we see that £khA.rt = (kr c ) ~ x l 2 . Since £khA_rt < 1 for A < r c , this implies that KH instabilities at any 
wavenumber will grow faster and are expected to play a larger role in the destruction of the clump than RT 
instabilities. 



The above timescales are for a hydrodynamic adiabatic shocked clump with an ideal gas equation of 
state. The expressions do not take into account the effects of additional physical processes which may be 
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present in shocked clump systems. Potential important processes include radiative cooling, self-gravity, heat 
conduction, and magnetic fields. By careful choice of physical parameters, we may include or exclude some 
or all of these. Magnetic fields, for example, are known to permeate much of the ISM. The role of magnetic 
fields therefore may be a significant one, and magnetic fields are an important aspect addressed in many 
studies (e.g. Jones et al.|1996 Shin et aL||2008 1 . Owing to the ubiquity of shocked clumps, however, we may 
consider the case of nonexistent or dynamically unimportant magnetic fields without loss of generality. 

The collapse of molecular cloud cores under their own gravity is understood to be one of the progenitors 
of star formation. While other factors may contribute, such as the passage of a supernova shock which 



compresses the cloud (Preibisch et al. 2002) or internal energetic jets which drive turbulence within the 



cloud Carroll et al. (20091, self-gravity is a universal component. The relevant timescale associated with the 



role of self-gravity is the free-fall timescale tg, 



ts 



(—) 
\32GpJ 

As discussed in § [3j our choice of scales is such that tg ^> t c 
ignored. 



1/2 



(4) 



and the effects of self-gravity may be safely 



As discussed by Orlando et al. (2005 1, thermal conduction may be ignored when the thermal conduction 
time scale t tc is cither much larger or smaller than the other physical time scales. Moreover, as described in 
KMC94 and others, one may postulate the presence of a magnetic field which inhibits thermal conduction 
while remaining dynamically unimportant. 



2.2. Radiative cooling 

A description of the role of radiative cooling in shocks may be found in texts such as |Dyson fc Williams] 



(1997); Zel'dovich & Raizer (20021. Shock waves convert kinetic energy directed perpendicular to the shock 



front into random thermal motions. Radiative cooling removes thermal energy from the system at a rate that 
is characterized by some cooling law. Analytic expressions for this cooling rate typically assume a power-law 
form of A(T) oc T a , where a = —1/2 may be assumed for typical ranges of temperatures (10 5 -10 7 K). From 



e.g. Dyson & Williams (19971, the cooling occurs over a time, 



3k B T ps v 3 s 

tcool = , oc — (5) 

where the subscript ps refers to postshock values and n a is the preshock (ambient) number density. We 
have made the usual assumption that T ps = p ps /(p ps kB) — {Pavf) / '(^Paks) in the strong-shock limit. The 
corresponding cooling length is 



Vpstcool 



!5 

n„ 



(6) 



where again we use the strong-shock relation v ps — 3w s /4. Since the pressure remains roughly constant over 
the cooling length, the ideal gas equation of state dictates that the density must increase. Hence, cooling 
shocks allow density contrasts greater than the factor of 4 expected from the strong shock limit of jump 
conditions with 7 = 5/3 ( Dyson fc Williams||1997 1. 
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It is often useful to consider a system in the isothermal (7=1) regime, which may be considered a 



limiting case of radiative cooling. From Dyson & Williams ( 1997 1, this regime occurs when L coo i — > and the 



temperature increase across the shock vanishes. In this regime v ps — v s instead of v ps = 3w s /4. Therefore, 
for a fixed postshock velocity (such as a specified boundary condition), the corresponding shock velocity will 
be reduced in the isothermal regime compared to the adiabatic regime. We therefore expect a decrease in 
the shock velocity as cooling strength increases. 

There are several visual cues to the cooling strength, one of the most obvious being a greatly reduced 
bow shock stand-off distance for a shocked clump. We may approximate this as follows. A strong shock 
impinging on a dense body will form a reverse-facing bow shock whose shape is determined by the Mach 
number of the flow and the shape of the solid body. In the adiabatic case with a rigid (\ = p c / Pa = 00) 
spherical body, the distance Aq between the bow shock and the body is ( Hirschei||2005| ) 



A . £l (7) 

where 

= p ± _ (7-1)^ + 2 
£ ~p 2 ( 7 +l)Af 1 2 (ti) 
with pi, p 2 the pre- and postshock density, respectively, and Mi the Mach number with respect to the 
preshock medium. For M >■ 1, e — > (7 — l)/(7 + 1). Thus, for 7 = 5/3, e — > 1/4 and the stand-off distance 
becomes A = l/3r c . 

What should happen to the bow shock stand-off distance when cooling is present? We approximate the 
effect of cooling by again considering 7 ~ 1. (This approximation works well in regions of compression; it 
breaks down in rarefaction regions — see Cunningham et al. (2009).) By inspection of Eq. [7] and [8] we see 



that in the 7 — > 1, M 3> 1 limit, Ao — * 0. We find therefore that as cooling strength increases, the bow 
shock stand-off distance Ao should decrease. This rough analysis is supported by numerical simulations of 
cooling shocked clumps: when cooling is present, the bow shock stand-off distance decreases compared to 
the adiabatic case. Thus strong cooling is characterized by a bow shock wraped tightly around the clump 
which is, itself, undergoing strong shock processing. 

We note that the situation is complicated somewhat when the rigid body is replaced by a solid body 
with finite density contrast x (i- e -j a clump). The bow shock reaches its steady-state stand-off distance in a 
time t os = Aq/w&s, where v^, s is the bow shock velocity. As soon as the clump is hit by the shock, however, the 
transmitted shock begins to accelerate clump material. The bow shock may not be accelerated downstream 
at the same rate as the post-transmitted shock material, since the stand-off distance is determined partly by 



the shape of the clump, which is changing. As the stand-off distance is greater for blunt objects (Hirschel 



2005), the axial compression of the shocked clump may push the bow shock further upstream from the clump, 
reducing its the bow shock's acceleration relative to the clump. Therefore, if t os < t cc , the bow shock may 
achieve a stand-off distance A os > Ao- (Conversely, if tf, s 3> t cc , then A os < Ao, and the bow shock may 
not achieve the steady state stand-off distance Ao given above.) 

In practice, in the adiabatic case of a shocked clump with finite x, tf, a < t cc , and A(, s reaches a steady 
state value and then remains stationary (in the preshock rest frame) for 1-2 t cc before being accelerated 
downstream. When cooling is present for the bow shock material, however, the reduced Vb s (leading to 
increased tb s ) results in a bow shock which forms a small distance from the clump, and remains tightly 
wrapped around the clump material over the course of the destruction of the clump. 



One therefore needs to define when the effects of cooling are important. In order for cooling to be 



- 7 - 



dynamically significant, the cooling length should at most be L coo i ~ r c . When L coo i < r c , cooling occurs 
rapidly and can be expected to be of greater importantance in the dynamics of the flow. As in |Fragile et al.| 
( [2004] ), we modify Eq. [5] for the case of a shocked clump using \ an d the shock jump conditions as, 

v 3 M 3 T 3/2 

tcooi = /3^7f— = 1.1 x 10 2 . ° s (9) 

where the subscript a refers to the ambient (preshock) material, and f3 = 6.61 x 1CP 11 cm~ 6 s 4 . 
Similarly we write a cooling length for the clump material as 

Lcooi — Vp S c t coo i — ~t CO ol 

Vx 

v 4 M 4 T 2 
= /Hr- = 1.3 x 10 6 ^^- cm (10) 

X n c X 3 n a 

where v pSlC is the postshock velocity in the clump. 

From a numerical standpoint, adequate resolution of the cooling length is important to ensure accurate 
capturing of the energy loss as the postshock fluid parcels cool. Inadequate resolution may lead to an under- 
estimation of the cooling rate or, cquivalcntly, an overestimation of the cooling length. In particular, the 
contact discontinuity between the upstream-facing clump bow shock and the downstream-facing transmitted 
shock must be adequately resolved in order to properly allow the growth of instabilities in the ablation 
region. This is discussed further in §§[5] and [6] 



Numerics 



3.1. The AstroBEAR code and the computational domain 



For the simulations we employed the AstroBEAR code using a 2D computational grid including the 
effects of axisymmetry (2.5D) as an operator-split source term. Employing a 2D grid with such a source 
term is widely employed in situations where variation in the polar angle <f> is expected to be minimal. 
Such is the case here with an initially planar shock propagating in z and a spherical clump located on the 
axis. We may therefore utilize a 2.5D simulation to approximate the results of a fully 3D simulation at 
reduced computational cost. However, to support our conclusions from the 2.5D data, we have undertaken 
a supplementary 3D simulation. This is discussed in § [4j 

The AstroBEAR code is a parallel AMR Eulerian hydrodynamics code with capabilities for MHD in two- 
and three-dimensions. There are several schemes of varying order available for the user; here the code solves 
the fluid equations using a second-order MUSCL scheme and treats cooling and cylindrical symmetry as 
Strang-split stiff source terms. AMR is employed to more accurately track regions of interest, determined by 



gradients in fluid variables. Details on AstroBEAR may be found in Cunningham et al. (2009 1. The inclusion 



of simple multiphysics is possible, tracking helium and hydrogen ionization; here hydrogen is assumed. The 



cooling curve of Dalgarno & McCray ( 1972 1 is employed; refinements on this method are being undertaken. 



Advected tracers are used to track the preshock, postshock, and clump material throughout the simulation. 

The computational domain has extent in z and r of 8r c and 2r c , respectively. Our computational domain 
is the same as that in Fragile et al. ( 2004[ ) and comparable to the 3D domain of Orlando et al. (2005). The 
physical parameters are listed in Table [T] The clump radius r c is 50 AU — however, this length scale is 
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not restrictive as our results scale to larger or smaller cloud sizes via appropriate choice of dimensionless 
numbers; sec below. The clump is located on the axis of symmetry at z = 2r c and has a density contrast 
with the ambient (preshock) medium of \ = Pel Pa = 10 2 - The preshock temperature is T a = 2000 K. A 
shock of mach number M = 50 was placed initially at z = 0.5r c , with postshock conditions imposed on 
the left (z = 0) boundary. The top and right boundaries have Oth-order extrapolating boundary conditions. 
The initially imposed shock crosses the domain in ~ l.lt cc . The time for all the clump material to blow 
off the grid is longer: the simulations all continue until all clump material has left the grid which occurs at 
approximately 3t cc . 



Quantity 


Ambient 


Clump 


Density 


10 2 cm" 3 


10 4 cra~ 3 






x= io 2 


Temperature 


2000 K 


20 K 


Sound speed 


5.2 km s _1 


0.52 km s" 1 


Shock velocity 


260 km s" 1 


26 km s _1 


Shock Mach number 


50 


50 


Radius 




50 AU 


tcooi (M=50) 


3.6 yr 


0.036 yr 


tcool 1 tec 


0.40 


0.004 






1.8 AU (0.04 r c ) 


ta 




5 x 10 5 yr 



Table 1: Physical parameters of the simulations for the ambient (preshock) and clump material. The shock 
velocities are for the external shock and transmitted (internal) shock, respectively. 

The resolutions for the 8 runs are given in Table [2] AMR was employed for the higher resolution runs, 
up to a maximum of 2 additional levels, where each subsequent level doubles the resolution. The "effective 
resolution" listed is the equivalent fixed-grid resolution for the highest refined level. Our lowest-resolution 
run had N=12 cells/r c . Each run thereafter doubled the resolution, to N=24, 48, 96, 192, 384, 768, and 1536 
cells/r c . We refer to each run according to these values as Rn- 



Run 


Effective Resolution 


Ax/L coo i 


Z\x' j L coo l .in spec. 


cells /r c 


cells in z x r 




Approx. 


12 


96 x 24 






24 


192 x 48 


1 


~ 1 


48 


384 x 96 


2 


1-2 


96 


768 x 192 


4 


3 


192 


1,536 x 384 


7 


5 


384 


3, 072 x 768 


15 


7 


768 


6,144 x 1,536 


29 


11 


1536 


12,288 x 3,072 


58 


14 



Table 2: Resolutions of the 8 simulations, from 24 to 1,536 cells per clump radius r c . The third column gives 



the number of cells per cooling length L coo i from Eq. 10 and the fourth column via visual inspection (see 
§ [B]) . At i?i2 the cooling length was unresolved. 

We note that the specific scales we have chosen are not critical for our results, with the exception of 
the preshock temperature, as the cooling rate is a function of temperature. The scales of importance are 
the dimensionless ratios such as the cooling length versus the clump radius L coo i/r c , the cooling time versus 
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the cloud-crushing time t coo i/t cc , and the density contrat x- We expect the results presented here to be 
demonstrative of any environment with similar dimensionless numbers. For comparison, for our choice of 
parameters, t cool /t cc = 0.004, identical to case E5 of |Fragile et al. (2004) which used a cloud of radius 
R c ~ lOOpc. 

The cooling length for the postshock clump material (see §[2]) is L coo i < 0.1r c , indicating that cooling is 
important in the clump. Since the preshock temperature is T a = 2000 K, cooling also will be important for 
shocks in the ambient material. As we will see, the choice of temperatures results in the strongest cooling 
occuring in the clump bow shock. 

We choose a sharp cloud boundary for our simulations, such that at radius r < r c the clump has 
density p c and temperature T c , and at radius r > r c these values are the preshock density and temperature. 
This yields a step-function profile. The choice of clump profile is a matter of some interest. Most studies 
have adopted a sharp boundary(e.g. KMC94, |Fragile et aT|2004| |Orlando et~al"1|2005| |Mellema et~aT|2002 



and others), while some recent studies have chosen to enclose a dense "core" within a "cloud" of larger 
radius which has some decreasing profile in density and corresponding transition in temperature to the 
preshock ambient medium (Poludnenko et al. 2002 Patnaude & Fesen 2005 Nakamura et al. 2006). As 
seen in Nakamura et al. ( |2006| ), the chief dynamic result of this is to soften the impact of the shock and 
increase the growth time of instabilities. We adopt a sharp boundary for conceptual simplicity. Such a 
boundary imposes a grid-based set of perturbations on the clump interface which in principle will influence 
the growth of Kelvin-Helmholtz (KH) and Rayleigh- Taylor (RT) instabilities. However, these perturbations 
have wavelengths A <C r c . Under the assumption of § [2] that the dominant modes are those with A ~ r c , 
these perturbations therefore should have minimal influence on how the solution varies by resolution. This 
assumption appears borne out by the results of our 3D simulation, discussed in § [4] which feature little or 
no grid-based artifacts. 



4. Results 

The same basic physics applies for all runs; see Fig. [T] A Mach 50 shock moving in the positive z 
direction with velocity v s impacts the clump located on the axis at z = 2r c with density contrast x = W0. 
This impact creates an internal (transmitted) shock which propagates through the clump at a slower velocity 
v s , c ~ v sX~ 1 ^ 2 - Both the internal and external shocks are radiatively cooling. Downstream (positive z) is 
rightward on the page; upstream is to the left. 

In a time t ~ t cc the external shock propagates off the edge of the computational domain (located at 
z = Sr c ). A reverse-facing bow shock forms which has a small stand-off distance from the edge of the clump 



(see § 2.2). This implies that the bow shock/contact discontinuity region is resolved by a correspondingly 
small number of cells. This important point is discussed further in §|5] Note as the internal shock propagates 
through the clump, material is stripped off the leading edge of the clump by KH instabilities and the bulk 
post-shock flow. We refer to this region (between the bowshock and the contact discontinuity) as the ablation 
front. RT bubbles propagate into the clump, while RT fingers are destroyed by the postshock flow. Finally, 
the downstream flow exhibits multiple interacting shocks. 

We now visually compare the eight different resolutions. Figure [2] gives synthetic Schlicrcn images of the 
simulations. Synthetic Schlieren images are a useful tool for investigating dynamic systems of shocks because 
they highlight the locations of gradients in density. Measured in cells per clump radius, the resolutions are 
12, 24, 48, 96, 192, 384, 768, and 1536. Each panel in the figure is a different resolution, and all panels are 
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Fig. 1. — Important features of the flow are given in a synthetic Schlieren image of run -R1536 at t ~ 0.5£ cc . The image 
has been reflected about the axis of symmetry, with the reflection showing the location of AMR refined regions. The 
bow shock wraps tightly around the clump from the strong cooling. The clump surface is ablated by its interaction 
with the postshock flow. A slip stream forms behind the clump. Transmitted shocks propagate internally through 
the clump. The external shock is susceptible to the cooling instability, and a conical reflected shock forms off of the 
axis of symmetry which reengages the flow. 



Fig. 2. — Synthetic Schlieren images of the 8 simulations of radiatively cooling shocked clumps in 2.5D are shown at 
time t = 1.44t cc . The Mach number of the shock is M — 50. Spatial units are in clump radii r c . A clump was located 
initially at z = 2r c and has a density contrast with the preshock ambient medium of % = 100. The computational 
domain in z and r is 8r c by 2r c , respectively - in the images the computational domain has been truncated on the 
left edge at z = r c . Each panel is labelled with the initial number of cells per clump radius (see § [3] and Table [2J. 
Visually, the simulations appear to show self-convergence up to a point (384 or 768 cells/r c ), after which they do not. 
The parameter regime is such that both the external and transmitted shocks are cooling (hence the extremely small 
clump bow shock stand-off distance, see §[2}. See §12] for more details. 
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at the same time in the simulation, t — 1.44t cc . (Note that in the image the left border of each panel is 
truncated at z = r c .) At this point in the simulation, the internal shock has completely crossed the clump. 

Are the simulations converging? In a typical demonstration of self-convergence, authors will show qual- 
itatively and quantitatively that as resolution increases, differences between successive resolutions decrease. 
We discuss quantitative apects the simulations and their convergence in § |5j Qualitatively, with increasing 
resolution there are visual cues of convergence, namely the presence of structures whose appearance remain 
basically the same, but "sharpen" with higher resolution. The ablation front remains at roughly the same 
location in z, while the width (extent in z) of the feature decreases with resolution. There is an extended 
feature starting at the rear of the clump (z = 3.0, r — 1.0) which extends downstream off the edge of the 
grid. This feature is a result of a conically shaped shock reflecting off the symmetry axis and reengaging 
material flowing around the clump. This results in a multiple-shock, axially-aligned structure. In cases 
through i?768, this structure is recognizable and is sharpened similar to the ablation front. At -R1536, this 
structure is destroyed apparently due to enhanced vorticity. Indeed, -R1536 appears to show the most dis- 
parate structures of any of the cases. This suggests that, while the simulations appear to converge up to a 
point, above a certain resolution key features of the simulations diverge as new physically relevent length 
scales are captured. 



2 4 6 80 2 4 6 8 




2 4 6 80 2 4 6 8 



Fig. 3. — Eight panels a-h show false-color representations of the logarithm of density for run -R1536 at eight times 
in the simulation, at f=-0.04, 0.4, 0.8, 1.2, 1.6, 2.0, 2.4, and 2.8 t cc , respectively. Both external and internal shocks 
are radiatively cooling. The clump bow shock is wrapped tightly around the clump. All of the clump material is in 
motion by the fourth panel, and RT and KH instabilities destroy the clump after a few t cc . See §|4]for further details. 

How important is this apparent disparity of -R1536? In other words, how does the evolution of case 
^?i536 compare with a case having resolution which in the adiabatic regime would be considered "sufficient"? 
We address this question by examining the time evolution of case -R1536, and comparing this evolution with 
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case i?i92, which is at a resolution that would be considered sufficient for an adiabatic simulation. Figure [3] 
follows the time evolution of the highest resolution run, -R1536, with a false color representation of logarithmic 
density. The eight panels a-h depict the simulation from t=-0.04 to 2.8 t cc . Panel a of Fig. [3] shows the initial 
conditions for the highest-resolution simulation. In panel 6, at t = 0.40i cc , the external shock has reached 
z = 4r c . The shock is not a sharp discontinuity owing to cooling. The internal shock has propagated ~ 2/3r c 
into the clump. An upstream-facing bow shock has formed around the clump, and is wrapped tightly around 
the clump surface. The surface of the bow shock at distances greater than r c from the cloud center is smooth. 
The ablation surface of the clump is being disrupted by RT and KH instabilities. The amplitude of the RT 
bubbles increases with radius, suggesting that there is a KH component which enhances the growth: further 
from the axis, the post-shock flow is increasingly oblique to the clump surface, implying that the relative 
velocity is higher and so from Eq. [2] the KH growth rate is increased. It is worth noting that there is a clear 
separation between the internal shock at z = 1.7 and the ablation front at z — 1.5. This is in contrast to 
other studies where the cooling rate is so high that the distance between these features essentially goes to 
zero (e.g. Mellema et al.|2002 Fragile et al.|2 004). This important point is discussed further in §Sj2]&|5] We 



also note the presence of a high-density peak at roughly z = 1.0, r = 0.8, whose evolution will be followed 
over the next few panels for comparison with i?ig2- 

In panel c at t — 0.8i cc , the external shock has reached z = 6.5 and continues to be disrupted. The 
internal shock has propagated about 75% of the way through the clump. The portion of the external shock 
seen wrapping around the clump in panel b has by this point reflected off the axis of symmetry, and the 
accompanying pressure increase has formed an upstream-propagating internal shock which has propagated 
into the rear clump by about 20%. The high-density peak noted in panel b has become separated from 
the clump. Though not readily apparent in this figure, its separation is caused by a fast-growing RT 
instability which forms at a radius just inside that of the peak, which then quickly propagates through 
the clump. The RT bubble is evident in the small downstream-facing shock emerging from the clump at 
roughly z = 3.0, r = 0.7. Overall, the downstream-facing internal shock front remains smooth in directions 
perpendicular to the flow, though its shape is distorted (especially at larger radii) by the instabilities and 
other transmitcd shocks compressing the clump. The front of the clump continues to be processed by the 
post-shock flow via the KH and RT instabilities. 

Panel d depicts t = l.2t cc . By this point the external shock has propagated off the grid. The region 
downstream of the clump exhibits KH instabilities, as well as the presence of numerous shocks. The shocks 
are a result of reflections off the axis as well as episodes of shock generation upstream near the clump. 
(The small shock noted in the previous panel is an example of one such episode.) There are also two large 
rarefaction regions downstream of the clump. The high-density peak noted in panel b, which separated 
from the clump in panel c, has further fragmented into two distinct vortical features in this panel, located 
at (z, r) = (3.2, 1.4) and (3.6, 1.7), respectively. The downstream-facing internal shock has propagated 
entirely through the clump and has collided with the upstream-propagating internal shocks. The shape of 
the clump ablation front is noticeably different than in earlier panels. The instabilities at the front have 
allowed a large portion of clump material to shear radially, leaving only material close to the axis behind. 
All of the clump material is in motion. 

Panels e through h, at times t — 1.6, 2.0, 2.4, & 2.8 t cc , respectively, exhibit the continued processing 
of the clump by the postshock flow. Clump material is episodically, rather than continually, stripped from 
the clump. The shape of the bow shock changes during these episodes as it remains tightly wrapped around 
the clump material. 

We now consider the evolution of case -R192, in comparison to that of i?i536- Figure [4] shows eight images 
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Fig. 4. — Eight panels show false-color representations of the logarithm of density for run R192 at eight times in the 
simulation, at t—-0.04, 0.4, 0.8, 1.2, 1.6, 2.0, 2.4, and 2.8 t cc , respectively. The figure may be compred with Fig. [3] 
While overall the evolution is similar there are notable differences, including the shape of the internal shock and the 
lack of vortices in the downstream region. 



- 15 - 



of case -R192 at the same times as shown in Figure [3] Overall, the behavior is similar. In particular, the 
axial location of the ablation front is consistent with that seen in Fig. [3] The high-density clump parcel 
noted in Fig. [3] is present and evolves similarly, appearing in panel c at z = 2.4, r = 0.8. However, in this 
case it remains a single vortex rather than multiple vortices. Also in panel c, the shape of the internal shock 
is different, as are the RT bubbles behind the clump ablation front. Here, the shock front remains smooth 
and unaffected by the RT bubbles up to the time when the shock crosses the clump, whereas in Fig. [3] the 
instabilities propagated up to the shock front, causing it to distort from planar. Overall, fewer vortices 
are noted, especially in the downstream region. We conclude that, visually, the behavior is unclear as to 
convergence, and deserves more quantitative investigation (§ 5). 

We undertook a single cooling 3D run for comparison with the 2.5D results, i?i92,3D, depicted in Fig. [5] 
Both panels in the figure are at the same time, (t = 0.8i cc ). The top panel depicts the logarithm of density 
in a volumetric plot. The simulation followed the ±x/+y/+z quadrant of the clump sphere, which has been 
reflected around to the -y/+z and -y/-z quadrants in the figure. The internal shock is evident in the image, as 
are the instabilities propagating into the clump. Note however that the 3D nature of the simulations allows 
us to capture non-axisymetric flow patterns, and the cloud surface is seen to be ablating in a nonuniform 
manner with respect to </>. Streams of material are evident flowing downstream from non-axisymetric clumps 
on the cloud surface. Information about the surface of the clump is propagated downstream as increased 
or decreased axial flow downstream. These non-axisymetric modes may have been seeded due to grid-based 
effects from the sharp cloud boundary. It is clear however that by this time in the evolution their effects 
on the cloud surface are minimal. If grid-based effects played a large role in the evolution one would expect 
corresponding features on the order of Ax. In fact, all visible perturbations at this time are appreciably 
larger, suggesting that the sharp cloud boundary does not play a significant role after providing an initial 
perturbation seed. 

The bottom panel depicts a "slice" through the computational domain in the y=0 plane of Rig2,3D 
compared directly to -R192 at the same time, which gives the same resolution but in 2.5D. A high degree of 
agreement is seen. The shapes and locations of the internal shock are nearly identical. The ablation fronts 
are very similar, though the 2.5D case has a more-developed RT bubble at z = 2.0, r = 0.6 than the 3D 
case. The rear of the clump, and the other transmitted shocks, also are quite similar. The shapes of the bow 
shock at radii larger than r c differ slightly, with the 2.5D case having a smoother rollover than the 3D case. 
Finally, the conical reflected shock is more prominent in the 2.5D case than in 3D. This good agreement 
suggests that the 2.5D simulations are reasonably reliable. 



5. Analysis 
5.1. Convergence of global quantities 



We wish to evaluate quantitatively how the simulations change with resolution; namely, whether they 
converge. We begin by examining the behavior of global quantities related to the clumps. A typical character- 
ization method of shocked clump simulations is to integrate certain physical quantities over the computational 
domain, and witness how these quantities vary with time (e.g. Nakamura et al.|2006 Pittard et al.|2009 Shin 
et al.|2008 1. KMC94 employed a two-fluid model, with one fluid being the clump and the other fluid the pre- 
and postshock medium, allowing them to explicitly track the clump fluid parcels. Nakamura et al. (20061 



employed a single-fluid model and used a density cutoff to prescribe "clump" material. Several authors track 
the quantities by including color (advected) tracers to differentiate the clump from the ambient environment 
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Fig. 5. — Top panel: A volumetric representation of logarithm of density for run Rig2,3D at t = 0.80t cc . The 
simulated region occupied one quadrant of space only; it has been reflected to two other quadrants in the image. 
The location of the transmitted shock is evident, as are the RT and KH instabilities at the clump ablation surface. 
Streams of material stripped off the clump surface are evident as well. Bottom panel: A comparison between runs 
-R192 and -Ri92,3D- A slice through the 3D domain of log density is shown in comparison to the 2.5D result at the 
same time, t = 0.80i cc , corresponding to panel c of Fig. [4] The results are in good agreement, differing slightly in 
the growth of RT bubbles into the cloud and the shape of the rollover region at the top of the clump. This suggests 
good reliability of the rest of the simulations presented. 
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(e.g. Pittard et al.| [2009 Shin ct al. 2008). Our advected tracers allows us to repeat this characterization. 



We define six global quantities as in Nakamura et al. (20061: clump mass, average axial velocity, radius 



in the cylindrical direction, radius in the axial direction, and velocity dispersions in the cylindrical and axial 
directions, respectively. 

The clump mass is defined as 

m c = [ P dV (11) 

Jc>c min 

where the integration is performed only over cells which have a fraction of clump material above a cutoff 
C > C m i n , < G < 1. For the results presented in §[5] we use C m i n — 0.1 — however we find the behavior of 
most quantities robust to the specific choice of C m i n with the obvious exception of m c . 

The radial and axial radii of the clump are defined as 

1/2 



» = I l(r 2 ) 



2 

2\-|l/2 



(12) 



c=[5«z 2 )-(z> 2 )]^ , (13) 
and the radial and axial velocity dispersions are similarly defined as 

SVr = (V 2 r ) 1/2 , (14) 

Sv z =((v 2 )~{v z r) 1/2 , (15) 

where for a quantity q n , 

(q n ) = m- 1 [ pq n dV . 

Jc>c min 

Finally, the average axial velocity is (v z ). 

To quantitatively address convergence, we use the standard definition of relative error, 

s(M,N)= IQn - Qm1 (16) 
Qn 

where Qn (Qm) is a global quantity defined as above for run Rn (Rm), and N > M. 
Convergence implies that, given a resolution Rn, 

e(N, 2N) < e(N/2, N) . (17) 

That is, if simulations with higher resolutions converge on the solution, the errors between higher and higher 
resolutions should decrease. This implies that a progression of relative errors, showing monotonic decrease, 
is necessary to determine convergence. Conversely, a lack of monotonic decrease in e(M, N) may suggest at 
least two possibilities. One possibility is that the simulations may be of such inadequate resolution that an 
increase in resolution from Rjy to R2N returns a solution functionally as poor as that of i?j\r- On the other 
hand, it may be simply that the simulations are not converging. Such would be the case, for example, in a 
highly nonlinear or fully turbulent simulation, where even at a given resolution, minor changes to the initial 
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conditions would yield differing evolution even if statistical measures were steady. A third possibility, that 
"plateaus" of convergence may exist, is discussed §|6] 

We may therefore ask how the relative error e(M,N) varies over our set of simulations, and whether 
through a progression of increasing resolution from Rm to R 2 m results in a monotonic decrease of e(M 7 2M). 
Table[3]lists the relative errors between all pairs of runs for one global quantity, clump mass, as defined above, 
at time t ~ t cc . At this point in the simulation, the external shock has completely passed the clump, and the 
transmitted shock has propagated roughly halfway through the clump. The columns in Table [3] represent 
runs at resolutions ./V and the rows represent runs at resolutions M . The leading diagonal, then, shows the 
e(M,2M) pairs of interest. Note that this is a representative time fairly early in the simulation, t ~ t cc ; at 



later times the values differ sec 5 5.2 



The values on the diagonal do not monotonically decrease, meaning that increases in resolution do not 
consistently reduce the relative error. Note moreover that a single value of e{M, N) is not a sufficient condi- 
tion to determine convergence. For example, e(12,24) = 0.06, but e(12, 1536) = 0.20. If only simulations at 
i?i2 and i?24 were performed, one might conclude based on e(12,24) that R\ 2 was converged, when in fact 
they both differ largely from i?i536- Both R 12 and R24 are so insufficiently resolved compared to -R1536 that 
e(12,24) means little on its own. (Of course, the derived value of L coo i would be an independent indicator 
that simulations at this resolution were not sufficiently resolved.) 

This illustrates why a progression of e(M,N) is necessary. In fact, for convergence to occur both 
e(N, 2N) and e(N, N max ) must be motonically decreasing and small (where by N max we mean the simulation 
of maximum resolution). A small e(N, N max ) ensures the goodness (accuracy) of the solution at resolution 
Rn- A small e(N, 2N) ensures that the goodness is real and not due to a selection effect: since the temporal 
evolution of the global quantities differs by resolution (as we will see below), it may be possible to choose 
an instant in time at which two resolutions seem to agree, misrepresent ative of the general behavior over 
time. It is unlikely, however, that at such a time both e(N, N max ) and e(N,2N) will be small, unless the 
simulation does in fact converge. 

In Table [3] we satisfy the criteria that e(N, N max ) decreases with increasing N. However, these data are 
from early in the simulation, t ~ t CCl and with the most well-behaved global quantity, m c . We do not see 
the same behavior with other quantities, and we see no convergence later in time, as we now discuss. 
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Table 3: Comparison of relative errors for each pair of simulations. The relative error is e(M, N) — \Qn — Qm\/Qn, 
where Qjv is a global quantity for run Rn and N > M. The quantity here is measured clump mass at t ~ t cc . 
The columns are runs Rn and the rows runs Rm. If convergence were being witnessed, the values along the leading 
diagonal ( e(M,2M) ) would be monotonically decreasing. In fact, a large increase in relative error is noted at 
e(192, 384), roughly the resolution at which the resolution typically is assumed to be sufficient. See §[5] for details. 
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5.2. The behavior of global quantities over time 
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Fig. 6. — Plots showing relative error e(iV, 1536) as denned in § 5.1 are presented at five times in the simulation, 
t = 0.20t cc , t = 0A0t cc , t — 0.80t cc , t = 1.2t cc , and t = 1.6t cc . Only very early in time does convergence appear; 
subsequent plots show little consistent indication of convergence. See § |5.2| 



It is common to show convergence results for e(N, N max ) for each global quantity at a representative 



time. In Nakamura et al. (20061, relative errors were investigated at one moment in time which, by visual 



inspection of Figures 2 and 3 of that work, was prior to the onset of any KH growth. Pittard et al. ( 2009J), 
in contrast, present convergence results at a time well after the onset of instabilities. Convergence results 
in terms of the relative error e(N, 1536) for the simulations are shown in Figure [6] The times depicted are 
t = 0.2i cc , t = 0At cc , t — 0.8t cc , t — 1.2i cc , t = 1.6t cc . These plots are representative; by inspection of 
the quantities over the course of the simulation we find it is only very early in the simulation, t < 0.5i cc 
that convergence seems apparent. After this point, there is no clear convergence in the measured quantities. 
Indeed, though not shown here, the behavior of the relative error varies widely over even very small timescales. 
The times from 0At cc to l.6t cc correspond to panels b-e in Figures [3] and [4] The last panel at 1.6i cc is close 
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in time to the comparison of Schlieren images in Fig. [2] 

In the first panel at t — 0.2t cc , all the quantities except the radial velocity dispersion Sv r show con- 
vergence in their steadily decreasing relative error with increasing resolution. The radial velocity dispersion 
shows an increase in relative error from i?gg to -R192, after which it decreases. In the second panel, at 
t = 0.40t cc , only clump mass m c shows monotonic decrease. The other quantities all show decreases in rela- 
tive error until some resolution, at which point the relative error increased. The bottom three panels neither 
show monotonic decrease in relative error nor relative errors as small as those seen in the first two panels. 
By 0.8t cc , the transmitted shock has not yet propagated entirely through the clump, though instabilities 
have been formed. 

The lack of observable convergence in the bottom panels of Fig. [^implies that either the global quantities 
are not a good measure of convergence after some disruption time, t > t^ isrupt , or that it is inappropriate 
to consider convergence after tdismpt due to nonlinearity. It is possible that convergence is observed only 
before any macroscopic growth of KH or RT instabilities, or before the onset of turbulence or large-scale 
vorticity. Since the KH and RT growth is resolution-specific, this suggests that the small-scale instabilities 
do have a global effect, namely in determining the time tdismpt- The growth rates also will be sensitive to 
the bow shock standoff distance: undcrrcsolving the region behind the bow shock may artificially damp the 
RT or KH instabilities. Thus, resolving this region is very important. The fact that cooling decreases the 
bow shock stand-off distance compared to adiabatic simulations implies that a higher resolution is necessary 
to achieve a well-converged simulation. 



5.3. Resolving the cooling length 

To get an idea of how well resolved the simulation is in the clump, we may consider the transmitted 
shock cooling length. We discuss the cooling length in § [2] where to derive an expression for L coo i, Eq. [6] we 



assume AcxT^ 1 / 2 . Using the expression for the cooling length in the clump, Eq.llOl we give the approximate 
number of cells per cooling length in column 3 of Table [l] 

In reality A is not a simple function of temperature. From a practical standpoint, one may estimate 
the resolution of cooling length by visually examining the number of cells of the cooling layer behind shocks 
in the simulation. Figure [7] shows line plots of the temperature on-axis for the simulations. Each curve has 
been displaced slightly to differentiate it. Simulations at R 12 and R24 did not sufficiently resolve the cooling 
length and are not shown. By inspection, the number of cells which captures the decrease from T rnax behind 
the shock to within 10% of T ps at each resolution is, i?4§: 1-2, Rqq: 3, -R192: 5, -R384: 7, i?76s: Hj -Ri536 : 
14. The apparent cooling length decreases with increasing resolution. These should be compared to the 
calculated values from Eq |10| which range from 1-58 for R48-R1536 ■ The primary cause of the discrepancy 
between the observed and predicted values is the deviation in the cooling curve from A cx T -1 / 2 , which was 
assumed in Eq. [10] 

The observed cooling length was lower than the predicted cooling length; this implies that a choice of 



resolution based solely on approximations like Eq. 10 may not be sufficient in a dynamic system. Without 



refinement criteria based on cooling, a simulation may not always adequately resolve the cooling layers 
behind shocks. Inadequate resolution of cooling layers increases the cooling time, effectively increasing the 
cooling length. Therefore, energy remains which should have been removed. On the other hand, successive 
shocks may drive the material into a more (or less) radiative regime, resulting in faster (slower) cooling. 
This confluence of effects implies that it may be difficult to recognize improperly resolved regions after the 
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Fig. 7. — Six plots depict the on-axis temperature profile (arb. units) at time t = 0A2t cc . Only the region near the 
transmitted shock is shown, and the curves have been displaced from one another for clarity. From left to right, the 
curves are for 7?48, Rge, R192, R3S4, R76S, and -R1536. As resolution increases, the temperature directly behind the 
shock increases, and the region of cooling decreases. See § |5.3| 
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fact. 

While one may choose a spatial resolution which resolves the cooling length initially, in a dynamic 
system this resolution may become insufficient. We propose a simple refinement criterion based on the 
cooling length behind shocks, L coo i, as 



NAx > L coo i = f3- 



(18) 



where v and n are the velocity and number density locally, j3 = 6.61 x 10 -11 cm -6 s 4 from Eq. [5J and Aa; 
the cell size. This therefore specifies a resolution of at least N cells per local cooling length. Based on our 
experience we suggest that N must be of order 10. Other refinement criteria related to coolling have been 
offered, for example, Niklaus et al. (20091 investigate refinement by overdensity, cooling time, and a term 
which depends on the enstrophy and rate of compression. Further development and comparison of refinement 
criteria will be illuminating. 



6. Discussion and Conclusion 



"Convergence" has a clear meaning when a problem is analytically tractable: a set of simulations are 
converging if an increase in resolution brings it closer to the analytic solution. Conversely, there are some 
cases where convergence is less well-defined, such as simulations of fully developed turbulence. For problems 
where there is no analytic solution, one still expects increasing resolution to improve the solution, because 
this for example brings the Reynolds number of the simulation closer to the actual value. (One may sidestep 



this by specifying a viscous scale, as in Pittard et al. (20091.) Simulations of problems with no analytic 



solution typically rely on a demonstration of self-convergence, which treats some high resolution simulation, 
believed by physical arguments to be more than sufficiently resolved, as "the solution." Lower resolution 
simulations in turn are compared to this "solution," and a resolution is chosen which balances accuracy with 
computational expense. 

Self-convergence for adiabatic shocked clumps shows a resolution of -Rioo~-R20o is sufficient to resolve the 



hydrodynamics of the problem (e.g. KMC94, Nakamura et al. (2006 1). It is easy to see that the inclusion 
of other physical processes may create different resolution requirements. For example, AMR simulations 



involving self-gravity may choose only to base refinement criteria off of the local Jeans length (e.g. Truelove 
et al.|1998 l. As seen in the present work, optically thin radiative cooling prescribes scales in the problem as 
well which relate to the local cooling length L coo i. 

In particular, if L coo i is highly underresolved, then wherever this is true, the simulation will be artificially 
moved into the isothermal regime. In Mellema et al. (2002) and |Fragile et al" (20041, this seems to be the 
case: v ps — v s in the clump, where L coo i/Ax <C 1 because L coo i/r c i ump < 0.1 and there are 100-200 cells per 
clump radius. Thus, there is no separation between the transmitted shock and the contact discontinuity. The 
vanishing of this separation results in artificially damping the growth of KH and RT instabilities, analogous 



to the results seen in the investigation of the NTSI in Klein & Woods (1998) and in contrast to the results 



presented here, which show an increase of KH and RT growth with resolution. It will be fruitful to observe 
the results of repeating the above simulations, but with the cooling length correctly resolved (requiring 



several thousand cells per clump radius in the case of Fragile et al. ( 2004 1) . The main difference between 
those results and ours is that they sufficiently resolved the hydrodynamics, and insufficiently resolved the 
cooling process, whereas in our case we progress from not resolving either process correctly to a regime in 
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which we believe we resolve both. 

Cooling enhances the growth of KH and RT instabilities (which follows from the fact that both depend 
on the thickness of the velocity/density shear layer, which decreases as L coo i decreases) and, as we've seen 
here, this has a global effect on the simulation: the breakup and mixing of the clump proceeds differently, 
and measured global quantities do not agree across resolutions. There is a shorter window (proportional to 
tKH and tux) in which the measurement of convergence via global measured quantities will succeed. Since, 
without explicit viscosity, the power spectrum of vorticity is a function of resolution, once vorticies have 
developed one would not expect a decent correspondence between different resolutions. One must therefore 
stipulate an alternate criterion to ensure correct resolution of physical processes throughout the simulation, 



such as the criterion based on the cooling length given in Eq. 18 Doing so will ensure the simulation 
remains correctly resolved and converging toward the correct solution, even when standard measures of such 
convergence do not apply. 

In the present suite of simulations, we seem to be approaching a converged solution. One might expect 
the solution to converge at -R3072 or i?6i44- However, there exists the possibility that, were the resolution to 
increase indefinitely, the solution would again diverge after crossing another physically-motivated resolution 
boundary. One could imagine several such "plateaus" of convergence, between which the solution appeared 
to diverge. It becomes all the more important, therefore, both to ensure the inclusion of correct physical 
processes, and to fully understand what resolution requirements are imposed by them. 
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